home *** CD-ROM | disk | FTP | other *** search
/ AOL File Library: 2,801 to 2,900 / aol-file-protocol-4400-2801-to-2900.zip / AOLDLs / C++ Files Library / FFT and Complex Numbers Library / FFTLibrary.sit / FFT library ƒ / FFT.cp < prev    next >
Text File  |  1995-01-15  |  2KB  |  87 lines
  1. //
  2. //    File: FFT.cp
  3. //    Copyright ⌐1994-1995 James Wilson
  4. //    All rights reserved.
  5. //    
  6. //    Revisions:
  7. //        8-2-94: Took out the constant multiplication part of arg and
  8. //        placed into halfArg. halfArg is multiplied by p to get arg.
  9. //        Removed the "f" suffix from 1.0 on line 34. Replaced the temp
  10. //        variable with the slightly more descriptive name "swap."
  11. //        7-29-94: Replaced all instances of floating-point division with
  12. //        multiplication by reciprocals, since multiplication is much
  13. //        faster on PowerPC.
  14. //    
  15. //    This is the implementation file for the Fast Fourier transform (FFT)
  16. //    algorithm. I have not made an extensive effort to comment every
  17. //    line of code, as that would lengthen this document by several more
  18. //    pages. I assume that the reader or user of this code has a decent
  19. //    understanding of both the DFT and the FFT.
  20.  
  21. #ifndef __FFT__
  22. #include "FFT.h"
  23. #endif
  24.  
  25. void FFT(const long n, const long nu, ComplexArray spectrumOut)
  26. {
  27.     long                    n2, nu1, i, k, p, x, tempIndex, denom;
  28.     double                    arg, reciprocalN, halfArg;
  29.     Complex                    exponent, swap;
  30.     
  31.     // Initalization process.
  32.     x = 1L;
  33.     n2 = n / 2L;
  34.     nu1 = nu - 1L;
  35.     k = 0L;
  36.     reciprocalN = 1.0 / n;
  37.     halfArg = twoPi * reciprocalN;
  38.     
  39.     // First step: calculate the Fourier transforms of the input values.
  40.     // The results are scrambled and are unscrambled at the end (see
  41.     // below).
  42.     while(x <= nu)
  43.     {
  44.         for(i = 1L; i <= n2; i++, k++)
  45.         {
  46.             denom = LongPower(2L, nu1);
  47.             p = BitReverse(k/denom, nu);
  48.             arg = halfArg * p;
  49.             
  50.             exponent.Real() = cos(arg);
  51.             exponent.Imag() = -sin(arg);
  52.             
  53.             tempIndex = k + n2;
  54.             
  55.             swap = exponent * spectrumOut[tempIndex];
  56.             spectrumOut[tempIndex] = spectrumOut[k] - swap;
  57.             spectrumOut[k] += swap;
  58.         }
  59.         
  60.         k += n2;
  61.         
  62.         if(k >= n)
  63.         {
  64.             x++;
  65.             n2 /= 2L;
  66.             nu1--;
  67.             k = 0L;
  68.         }
  69.     }
  70.     
  71.     // Last step: unscramble the results. This is done by swapping the
  72.     // values in spectrumOut[k] with spectrumOut[BitReverse(k, nu)].
  73.     while(k < n)
  74.     {
  75.         i = BitReverse(k, nu);
  76.         
  77.         if(i > k)
  78.         {
  79.             swap = spectrumOut[k];
  80.             spectrumOut[k] = spectrumOut[i];
  81.             spectrumOut[i] = swap;
  82.         }
  83.         
  84.         k++;
  85.     }
  86. }
  87.